Comparative study on free vibration analysis of rotating bi-directional functionally graded beams using multiple beam theories with uncertainty considerations

The present study investigates the free vibration behavior of rotating beams made of functionally graded materials (FGMs) with a tapered geometry. The material properties of the beams are characterized by an exponential distribution model. The stiffness and mass matrices of the beams are derived using the principle of virtual energy. These matrices are then evaluated using three different beam theories: Bernoulli–Euler (BE) or Classical Beam Theory (CBT), Timoshenko (T) or First-order Shear Deformation Theory (FSDT), and Reddy (R) or Third-order Shear Deformation Theory (TSDT). Additionally, the study incorporates uncertainties in the model parameters, including rotational velocity, beam material properties, and material distribution. The mean-centered second-order perturbation method is employed to account for the randomness of these properties. To ensure the robustness and accuracy of the probabilistic framework, numerical examples are presented, and the results are compared with those obtained through the Monte Carlo simulation technique. The investigation explores the impact of critical parameters, including material distribution, taper ratios, aspect ratio, hub radius, and rotational speed, on the natural frequencies of the beams is explored within the scope of this investigation. The outcomes are compared not only with previously published research findings but also with the results of 3-Dimensional Finite Element (3D-FE) simulations conducted using ANSYS to validate the model’s effectiveness. The comparisons demonstrate a strong agreement across all evaluations. Specifically, it is observed that for thick beams, the results obtained from FSDT and TSDT exhibit a greater agreement with the 3D-FE simulations compared to CBT. It is shown that the coefficient of variation (C.O.V.) of first mode eigenvalue of TSDT, FSDT and CBT are approximately identical for random rotational velocity and discernible deviations are noted in CBT compared to FSDT and TSDT in the case of random material properties. The findings suggest that TSDT outperforms FSDT by eliminating the need for a shear correction coefficient, thereby establishing its superiority in accurately predicting the natural frequencies of rotating, tapered beams composed of FGMs.


C h
Thickness taper ratio C.O.V Coefficient of variation {d s } Displacement vector E Young's modulus (Pa) E [omega]  Mean approximated natural frequencies g x Gradient index through the longitudinal direction g z Gradient index through the thickness direction h(X) Beam thickness or height as a function of X (m) h 0 Beam thickness at the beam root (m) h L Beam thickness at the beam free end (m) I 0 Beam second area moment of inertia at the beam root (m Axial displacement on the neutral axis (i.e., z = 0) U 1 , U 2 Axial displacements at nodes 1 and 2 respectively v i

Random variables vi
Mean value of the random variable V Volume (m 3 ) Var [ω]  First-order approximated variance of natural frequencies w(x, z, t) Lateral or flapping displacement field w 0 Lateral displacement on the neutral axis (i.e., z = 0) Cross-section rotations of FSDT at nodes 1 and 2 respectively ψ Cross-section slope at z = 0 in TSDT 1 , 2 Line slopes at z = 0 of TSDT at nodes 1 and 2 respectively ω Natural frequency (rad/s) ω (0) Zero-order term of natural frequency ω (1) ,i First-order term of natural frequency ω (2)  ,ii Second-order term of natural frequency ω Dimensionless natural frequency(ω ρ 0 A 0 L 4 /E 0 I 0 ) Rotating velocity (rad/s) ¯ Dimensionless rotating velocity(� ρ 0 A 0 L 4 /E 0 I 0 ) δ Virtual operator Mode shape vector Rotating beams can be used to simulate structures such as industrial fans, helicopter and propeller blades, wind and steam turbines, robot manipulators, and spinning space structures 1 .Understanding the dynamical behavior of rotating beams is critical during the early design stages to avoid resonance within the operational speed range and to evaluate rotating beam performance.These structural components are made of various material classes to satisfy various engineering design specifications.Because of their superior properties, rotating beams manufactured from composite materials have been widely applied in a range of industrial applications during the last few decades.When compared to traditional engineering materials, the advantages of composite materials, together with their ability to customize their designs to specific uses, have given them a competitive advantage.Aksencer and Aydogdu 2 investigated the free vibration of a rotating laminated composite beam with attached mass using the Ritz method.The authors used different beam theories (CBT, FSDT, and TSDT) in the formulation and considered the cross-ply lamination configuration.In addition, attached mass to beam mass ratio, position of the mass and aspect ratio were examined.Based on Timoshenko beam theory and nonlinear strain-displacement, Khosravi et al. 3 discussed the thermal vibration of rotating composite beams.The beams were reinforced by Carbon Nano-Tubes (CNTs) according to a uniform and two symmetric gradient distributions.The authors used the differential transform method to obtain the natural frequencies for different parameters such as hub radius, rotating speed, aspect ratio, temperature, and boundary conditions.Xu et al. 4 used the TSDT with the Ritz technique to investigate the first resonance frequency of rotating nanocomposite beams reinforced with carbon nanotubes.Mohammadi et al. 5 used the Differential Quadrature Method (DQM) to solve the governing equations for a rotating beam composed of multilayer piezoelectric nanobeams.The authors employed surface elasticity theory in conjunction with nonlocal continuum theory for the Timoshenko beam to derive the equations of motion.Furthermore, they investigated the effects of four different boundary conditions, thermal stress, external voltage, and angular velocity on the natural frequency.
A class of composites called Functionally Graded Materials (FGMs) has attracted significant attention in several modern engineering applications.The required variational features of the material combinations are increased by the application of FGMs to increase functional performance.Because of this continual fluctuation, FGMs have a continuous stress distribution and prevent stress concentrations.Another exceptional quality of FGMs is their capacity to endure high temperatures while preserving structural integrity.Akbaş 6 used Lagrange's equations to derive the equations of motions and then the Finite Element Method (FEM) to obtain the thermal natural frequencies for axially FGM Bernoulli-Euler beam.The material properties were changed by a power-law and were temperature dependent.To examine the nonlinear free thermal vibration of pre/post buckled rotating FGM beams, Arvin et al. 7 employed Bernoulli-Euler theory and a nonlinear strain displacement relationship to propose some new algorithms in conjunction with the nonlinear finite element method.Van Dang 8 applied Bernoulli-Euler theory and the FEM to investigate the static bending of functionally graded porous rotating beams.The beam was impacted by lateral and axial compressive force and embedded in an elastic foundation with two parameters.The material properties vary in the thickness direction according to the power law.Binh et al. 9 obtained and solved the equations of motion for a rotating Timoshenko beam formed of functionally graded porous material reinforced by graphene platelets using the Chebyshev-Ritz method.The material characteristics change according to two different types of porosity distributions through thickness and two different graphene platelet dispersion patterns.The authors investigated the influence of rotational speed, hub radius, porosity, and weight fraction on the natural frequency.Dang et al. 10 considered the Coriolis and centrifugal forces in deriving the equations of motion by using Hamilton's principle with Love's shell theory for a cylindrical shell made of Functionally Graded Porous (FGP) material.The porosity was changed through the thickness according to three different porosity distributions.The authors used Galerkin's method to obtain the natural frequencies for different boundary conditions.Also, they studied the effect of porosity distribution, rotational speed, Coriolis acceleration, and geometric parameters.
Unfortunately, the 1D-FGMs are ineffective at meeting the technical specifications for shuttles and aerospace craft, such as the stress distributions and temperature in various directions 11 .Utilizing material properties that change in desirable directions, such as two-directional FGM, can address this limitation 12 .Fang et al. 13 investigated the time response and coupled axial, flap-wise, and chord-wise vibration of a rotating BFGM cantilever beam.The beam material was gradually changed according to a power-law though the width and thickness.Lagrange's equation and the Ritz method were used to derive the dynamic equations and then solved by using the state space method for different material gradients and rotating speeds.
Rotating beams can be classified based on their geometric properties as either uniform or tapered.The later is often preferred due to its ability to provide an optimal distribution of weight and strength, which is particularly useful in meeting specific structural and functional requirements 14 .Banerjee et al. 15 used Hamilton's principle to derive the equations of motion for the flap vibration of a rotating double tapered Bernoulli-Euler beam.The author used the Wittrick-Williams algorithm to solve the resulting dynamic stiffness matrix for different taper ratios, rotational speeds, and hub radii.Lagrange's form with the FEM were used to develop mass, elastic, and centrifugal stiffness matrices for a rotating tapered Bernoulli-Euler beam by Bazoune 14 .The author examined the tapering effect in two planes, hub radius, and rotational speed.Chen et al. 16 examined the accuracy and efficiency of the variational iteration method for the free vibration analysis of a rotating Timoshenko beam.The beam was linearly tapered though the width and height.Adair and Jaeger 17 reformulated the fourth-order differential equation as a first order matrix and used the power series method to obtain the natural frequencies of a rotating taper Bernoulli-Euler beam.The authors studied both cone and wedge cantilever beams for different taper ratios.Nourifar et al. 18 compared the differential transform method and the finite element method for the vibration of a rotating cylindrical tapered Bernoulli-Euler beam.The effect of rotating speed and taper ratio were examined on the natural frequencies.An improved transfer matrix method was developed by Lee and Lee 19 to obtain the bending natural frequencies for a tapered rotating Bernoulli-Euler beam.The Frobenius method for a power series was used to solve the equations of motion.The authors studied the effect of centrifugal axial force, taper ratio, and hub radius on the natural frequencies.Wang and Li 20 used the differential quadrature method to solve the differential equations obtained by Hamilton's principle for the lateral vibration of a tapered rotating hollow www.nature.com/scientificreports/beam.The beam has constant inner cross-section radius and tapered outer radius.The authors examined the effect of rotating speed, hub radius, aspect ratio, taper ratio, and inner radius.Some literature has explored the use of both taper and composite materials in rotating beams.Piovan and Sampaio 21 used the variational principle to develop a nonlinear model for a rotating FGM tapered Timoshenko beam.The FEM was used to obtain the natural frequencies for different material distribution, aspect ratios, and speeds.Zarrinzadeh et al. 22 conducted an in-depth investigation of the vibration characteristics of a rotating tapered axially functionally graded material (AFGM) Bernoulli-Euler beam using the finite element method.Their study encompassed a systematic exploration of various influential parameters, including material properties, taper ratio, rotational speed, hub radius, boundary conditions, and the presence of a tip mass.For both Bernoulli-Euler and Timoshenko rotating tapered FGM beams, Hajheidaria et al. 23 investigated the lead-lag, flap, and flap-lag vibration.The metal/ceramic-based FGM beam properties changed though the thickness according to a power-law in a symmetric structure.The authors used the finite element method with 4 DOF and 8 DOF element models.Also, they studied the effect of volume fraction, rotational speed, hub radius, and taper ratio.Kumar et al. 24 used the differential transform method to estimate the flap wise natural frequencies of tapered FGM beams.The material properties changed laterally from the middle to the outer surface symmetrically and were estimated using Mori Tanaka methods.The authors discussed the effect of rotational speed, hub radius, taper ratio, and gradient index on the frequency.
The modified variational method and multidomain mixed approximations were used by Tian et al. 25 to investigate the vibration analysis of a double-tapered rotating FGM beam.The beam material including porosities was distributed based on the modified rule of mixtures.In this model, the Coriolis and nonlinear effects were considered for bending-stretching, twist-stretching, and bending twist vibration modes.The authors investigated the material, rotation speed, and various geometric effects.Bhattacharya and Das 26 considered the non-linear geometry, Coriolis acceleration, spin-softening, and thermal environment to study the free vibration of rotating micro-beams.In this study, Timoshenko theory, along with modified couple stress theory, investigated a double taper BFGM rotating beam.Also, the authors examined the effect of FGM composition, size-dependence, taper ratio, aspect ratio, hub radius, and temperature.To obtain the natural frequencies and mode shapes of a rotating BFGM for a tapered cantilever beam, Zhou et al. 27 used the Rayleigh-Ritz method.The equations of motion for time-dependent rotating velocity with periodic coefficients were derived by using Hamilton's principle and the Galerkin method.Bolotin's method with a higher-order approximation is used to solve the dynamic instability caused by periodic rotational velocity.The authors examined the effect of hub radius, rotational speed, FGM index, dynamic amplitude factor, and taper ratio on dynamic instability and natural frequency.Özdemir 28 investigated the free vibration and buckling behavior of rotating beams.Considerations included linearly tapered beams and axially functionally graded materials using a simple power law.Bernoulli-Euler and Timoshenko beam theories were applied using the Finite Element Method.The author investigated many parameters including the hub radius, rotational speed, power law index, aspect ratio, taper ratio, and other boundary conditions.
Various beam theories have been employed to investigate the vibration characteristics of rotating beams.The classical beam theory, also known as Bernoulli-Euler theory (CBT), represents the oldest and most fundamental approach.CBT assumes that the cross-section of the beam remains planar and perpendicular to the beam axis after deformation.Due to its simplicity and suitability for thin beams where transverse shear deformation is less significant, CBT continues to be widely utilized [6][7][8]14,15,[17][18][19][20][22][23][24]27,29 .
Third-order shear deformation theory (TSDT), also known as Reddy beam theory, goes a step further, and accounts for the fact that the cross-section will no longer remains straight or perpendicular to the beam axis after deformation.In TSDT, the transverse shear strain and stress are assumed to have a parabolic distribution with respect to the thickness coordinate 4,8,40 .
In the literature, there is a paucity of studies that compare different beam theories for rotating beams.However, Hajheidaria et al. 23 and Özdemir 28 discussed and compared Bernoulli-Euler and Timoshenko beam theories.Furthermore, Aksencer and Aydogdu 2 conducted a comparison of the Reddy, Timoshenko, and Bernoulli-Euler beam theories for rotating beams.These studies provide important insights into the performance of different beam theories and can assist in the design and analysis of rotating beam structures.Nonetheless, further research is needed to fully understand the behavior of rotating beams and determine the most appropriate beam theory for specific applications.
The prior investigations have primarily focused on the vibration of rotating beams characterized by deterministic properties.Nevertheless, it is crucial to acknowledge that real-world structures and mechanical systems inherently possess random properties.These uncertainties have a significant impact on both performance and structural reliability.Within the realm of rotating structural systems, these uncertainties arise from various sources, including variations in loads and material properties.To design highly reliable rotating beam structures, it is essential to comprehensively examine the collective effects of uncertainties in material and sectional properties, geometric parameters, and angular velocity on the stochastic response of rotating beams.Furthermore, conducting sensitivity analyses is crucial as it allows us to identify the design parameters that have the most significant influence on the computed statistical characteristics of the structural responses.Several numerical studies have been conducted to explore uncertainties in the dynamics of rotating beams [41][42][43][44] .However, based on the existing literature and the authors' knowledge, no study has explored the uncertainty associated with material distribution in functionally graded material (FGM) rotating beams and its impact on their natural frequencies.

Mathematical model formulation
A bi-directional functionally graded material beam of total length L along the axial axis X and a double tapered cross-sectional area A(X) is shown in Fig. 1. b(X) is the width and parallel to the Y axis and h(X) is the thickness and parallel to the Z axis.The beam is attached to a rigid hub of radius R. The beam rotates in the X − Y plane with a constant angular speed in rad/s about the Z axis.Also, the beam is divided into N e elements with equal length ℓ and has a local coordinate system x, y, z , where X, Y , Z is the global coordinate system.L i denotes the offset of the ith element from the Z-axis as follows:

Geometry and material properties
The double tapered beam geometrical dimensions and cross-sectional area and moment of inertia is given as a function of x by (1) Schematic drawing of a double tapered rotating BFGM beam.
where x + L i = X and h 0 , b 0 , A 0 and I yy 0 are the thickness, width, cross-sectional area and second area moment of inertia at the beam root, respectively.h L and b L denote the beam thickness and width at the free end.Also, 0 ≤ C b = 1 − b L b 0 < 1 is the width taper ratio and 0 ≤ C h = 1 − h L h 0 < 1 is the thickness taper ratio.It is noted that C b = C h = 0 for a uniform beam cross section and for C b = C h the beam has different taper values in the width and thickness directions.
In this work, the beam material properties vary continuously along the beam thickness or the longitudinal direction or both according to exponential rule of mixtures 12 .Thus where P(x, z) is the effective material properties (Young's modulus E, density ρ and shear modulus G) and P 0 is the material property at the reference position 0, − h 0 2 .g x and g z are the gradient indexes through the longitu- dinal and thickness direction, respectively.The beam material is homogeneous for g x = g z = 0.

Displacement and strain fields
In the current work, the coupled axial and flap-wise transverse vibration is considered; hence the lag-wise or twisting vibration is not considered.The axial displacement u and transverse displacement w of any point on the beam according to CBT, FSDT, and TSDT are given by Eqs. ( 4)-( 6), respectively and shown in Fig. 2b-d respectively.Figure 2a represents the undeformed cross-section for reference.
The superscripts ( ) C , ( ) F , and ( ) T are used to represent CBT, FSDT, and TSDT, respectively.u 0 and w 0 are the axial and transverse deflection on the neutral axis (i.e., z = 0 ) and t is time.∂w 0 ∂x denotes the pure bending cross-section slope in CBT, φ is the cross-section rotation in FSDT and ψ is the slope of the deformed line at z = 0 in TSDT as shown in Fig. 2b-d respectively.
Equations ( 4)-( 6) can be rewritten in terms of the displacement vector d s as: (3) The cross-section deformation for CBT, FSDT, and TSDT theories 45 .
Vol where the superscript prime ( ) ′ denotes the derivative with respect to x and the superscript ⊺ denotes the trans- pose.For linear and small deformations, the non-zero strain fields can be represented by Eqs. ( 10)-( 12) for CBT, FSDT, and TSDT, respectively.
or in vector form as

Stress-strain constitutive equations
The stress-strain constitutive equations according to Hooke's law for linear and small deformations can be considered for the FGM beam for CBT, FSDT, and TSDT beam theory as, respectively: where κ s denotes the shear correction factor for FSDT.

Virtual energy expressions
The virtual potential energy expression due to the stress field of a single beam element can be obtained in the form where the subscript s indicates stress field, the superscript ⊺ denotes the transpose and V is the volume.Sub- stituting Eqs. ( 13)-( 18) into Eq.( 19) gives the virtual strain energy for CBT, FSDT, and TSDT, respectively, as The virtual potential energy expression due to the axial centrifugal force of a single beam element is given by where the subscript cf denotes centrifugal, and F cf is the centrifugal force that can be obtained as The virtual kinetic energy expression of a single beam element is where the superscript dot ( ) indicates the time derivative.The specific virtual kinetic energy for CBT, FSDT, and TSDT can be obtained by substituting Eqs. ( 7)-( 9) into the general form of virtual kinetic energy Eq. ( 25) as: where

Finite element modeling
The condition of equilibrium of a dynamical structure for free vibration based on the principle of virtual energy is To accurately analyze the vibration behavior of rotating double taper beams made of BFGM, FEM techniques are employed using three distinct theories: CBT, FSDT, and TSDT.Each theory requires a specific number of degrees of freedom (DOFs) to adequately capture the axial and transverse displacements, as well as the rotational motion of the beam.In the CBT approach, a two-node, six DOFs element is utilized to capture the desired displacement and rotation components.These DOFs comprise two DOFs for axial displacement, two for transverse displacement, and an additional two for rotational bending motion, as illustrated in Fig. 3a.Axial displacement is estimated using Lagrange linear shape functions, while the transverse displacement is estimated using Hermitian shape functions as given in Eq. (30).
To achieve higher precision in the analysis, the FSDT method is employed, utilizing a two-node, ten DOFs element to accurately capture axial and transverse displacements, as well as rotation.The ten DOFs element for Timoshenko beam analysis offers advantages over the six-DOF element, including a superior convergence rate 23 .The nodal displacement vector in this model consists of two degrees of freedom for axial displacement, two for transverse displacement, two for rotational motion, and their corresponding derivatives as shown in Fig. 3b.Axial displacement is estimated using Lagrange linear shape functions, while transverse displacement and rotation are estimated using Hermitian shape functions based on the nodal displacements as given in Eq. (31).
For enhanced accuracy in the analysis, TSDT theory is utilized, employing a two-node, eight DOFs element.These eight DOFs encompass two degrees of freedom for axial displacement, two for transverse displacement, two for the derivative of transverse displacement, and two for rotation as shown in Fig. 3c.The estimation of axial displacement and rotation relies on Lagrange linear shape functions, while transverse displacement is approximated using Hermitian shape functions, both formulated in terms of nodal displacements as given in Eq. ( 32).
where {q C } , {q F } , and {q T } represent the displacement, {q C e } , {q F e } , and {q T e } denote the nodal displacement, and [N C ] , [N F ] , and [N T ] represent the matrix of shape functions for CBT, FSDT, and TSDT, respectively and given in Appendix A.
By substituting Eqs. ( 30)- (32) into Eqs.( 20)-( 23) and Eqs. ( 26)-( 28), summing the energies for all elements, and then applying Eq. ( 29), we obtain the following equation of motion: where [M] is the mass matrix and [K] is the stiffness matrix of the beam.The free vibration analysis is applied to Eq. (33) to calculate the natural frequencies by solving the eigenvalue problem Eq. (34), where is the mode shape.

Finite element modeling using ANSYS
To validate the obtained results from the developed mathematical models, 3-Dimensional Finite Element (3D-FE) simulations were conducted using ANSYS Workbench.A model of a tapered beam was created using the ANSYS Design Modeler, see Fig. 4. The simulations were performed for both isotropic and functionally graded material beam models.The geometrical parameters under investigation such as taper ratios (C b and C h ) , aspect ratio L h 0 , and hub radius (R) were integrated into the model as variables to easily control them through the Workbench's parametric table.In addition to these parameters, other variables such as rotational speed ¯ and index parameters for the material distribution g x and g z were also included to assess their impact on the natural frequencies of the beams.
Quadratic hexahedral SOLID185 elements were used to discretize the tapered beam geometry.Mesh controls were utilized to manipulate the mesh size of the beam edges, ensuring that only hexahedral elements were present in the mesh.The implementation of functionally graded material properties in the model involved incorporating APDL commands within the modal module of Workbench.The material properties of the BFGM beam are determined at the coordinates of the centroid of each element, as described in Eq. ( 3).Thus, in order to achieve a uniform distribution of material properties, discretizing the beam into regular hexahedral elements is necessary.A typical mesh is shown in Fig. 4. The material properties were then assigned into the respective element.

Uncertainties and stochastic modal analysis
As commonly recognized, system uncertainties encompassing material properties (E, G, and ρ ) as well as mate- rial distribution ( g x and g z ) in a rotating beam are subject to fluctuations in proximity to their designated values throughout processes such as measurement, structural element manufacturing, and structure assembly.Consequently, it is imperative to consider system parameters as stochastic rather than deterministic.Furthermore, the dynamic behavior of a rotating beam differs from that of a non-rotating beam due to the additional influence of centrifugal forces.The angular velocity assumes a pivotal role in these phenomena and frequently exhibits variations in the vicinity of its operational speed.Thus, the angular velocity is inherently characterized by randomness.Henceforth, we will employ the notation v i ( i = 1, 2, . . ., 6 ) to represent the individual baseline random variables, with v 1 denoting , v i ( i = 2, 3, 4 ) representing the material properties E, G, and ρ , and v i ( i = 5, 6 ) signifying the material gradient indices g x and g z , respectively.These uncertainties have the potential to intro- duce variations in the components of the matrix [K] as outlined in Eq. (34).Given the uncertainty associated with the matrix [K], the natural frequencies themselves become stochastic variables.To ascertain the statistical characteristics of these natural frequencies, one can deduce them from the statistical properties of the baseline random parameters, employing a mean-centered second-order perturbation methodology.
The stochastic modal analysis of BFGM tapered rotating beams utilizes the mean-centered second moment method 41 .This perturbation approach is founded on the principle of expanding the random response around the mean values of the baseline random variables while retaining terms up to the second order.Within the framework of the mean-centered second-order method, the second-order approximate mean and the first-order approximated variance of natural frequencies can be mathematically expressed as follows: Here ω (0) represents the zero-order term of the natural frequency, equivalent to the deterministic natural frequency.ω (1) ,i signifies the first-order term of the natural frequency concerning the random variable v i , and ω (2) ,ii corresponds to the second-order term of the natural frequency with respect to the random variable v i .S ω and S v i are the standard deviation of ω and v i , respectively.Hence the first and second-order term of natural frequency according to the finite difference approximation: where v i = v i − vi , with vi denoting the mean value of the random variable v i .

Results and discussion
The axial and flap-wise bending vibration analysis of a rotating BFGM double-tapered cantilever beam is investigated for CBT, FSDT, and TSDT.The beam material varies according to the exponential law of distribution as given in Eq. (3).Several parameters such as the index parameters for the material distribution (g x andg z ) , rotating speed (�) , taper ratios (C b andC h ) , hub radius (R), and aspect ratio L h 0 are discussed.This section is divided into two subsections.The first subsection gives the convergence and validation of the present model with 3D-FE where ω is the dimensionless natural frequency, ¯ is the dimensionless rotating speed, and R is the dimension- less hub radius.
To investigate the statistical characteristics of natural frequencies in BFGM rotating beams with varying baseline random variables, we consider the following scenarios: Case 1: Only the rotational velocity ( ) is treated as a random variable.Case 2: Only the material properties (E, G, and ρ ) are considered as random variables.Case 3: Only the material distribution indices ( g x and g z ) are subject to randomness.
In order to obtain comparative results, we also employ the Monte Carlo (MC) method, generating 500 samples for analysis for case 1 only.

Convergence and validation
The current model includes many parameters that should be taken into consideration to ensure a through verification procedure.These various parameters are the gradient index parameters for the material distribution (g x and g z ) , taper ratios (C b andC h ) , rotating speed ( �) , hub radius ( R) , and aspect ratio L h 0 .

Convergence rate
The convergence rates of the first five dimensionless natural frequencies are presented in Table 1 for CBT, FSDT, and TSDT and compared with 3D-FE results, where the subscripts f and ax denote flapping and axial modes respectively.The following data is used: The analysis of the data presented in Table 1 indicates a rapid convergence of results as the num- ber of elements (N e ) increases for the CBT, FSDT, and TSDT beam theories.Notably, when N e reaches 100, the accuracy of the obtained results becomes satisfactory for practical purposes, and hence this number of elements is employed for conducting free vibration analysis of rotating BFGM double-tapered beams.Remarkably, the outcomes derived from TSDT, FSDT, and 3D-FE analysis demonstrate a noteworthy level of agreement, particularly for thick beams at higher modes, when compared to the classical beam theory (CBT).This disparity can be attributed to the inherent limitations of CBT, which neglects the significant influence of shear deformation effects that become increasingly prominent in the higher modes of thick beam structures.

Effect of gradient index and aspect ratio
The provided data in Table 2 presents the first three dimensionless natural frequencies of a non-rotating ¯ = 0 BFGM uniform beam (C b = C h = 0) .The beam is analyzed using three different theories: TSDT, FSDT, and CBT.The frequencies obtained from these theories are compared with the reference values calculated using the dimensionless parameter = ω √ 12 from 12 using TSDT and 3D-FE results.The beam dimensions and material properties used in the analysis are (b 0 = 0.5 m, h 0 = 1 m, ρ 0 = 7850 kg/m 3 , E 0 = 210 GPa, ν = 0.3, and R = 0) .Two different aspect ratios, L h 0 = 20 and 5, are considered.It can be observed from Table 2 that TSDT, FSDT, and 3D-FE demonstrate a commendable level of agreement with the reference results, regardless of whether the beams are thick or thin.This indicates the robustness and reliability of these methods in accurately predicting the natural frequencies.While CBT shows good agreement with the reference results for the 1st mode of a thin beam L h 0 = 20 , its performance deviates from the TSDT, FSDT, 3D-FE, and 12 results for other modes and thicker beams.However, it is worth noting that for the third mode of L h 0 = 5 , which corresponds to the first axial mode, all three theories (TSDT, FSDT, and CBT) exhibit satisfactory agreement with both the 3D-FE and reference results.
Furthermore, the results show that as the gradient index g x = g z increases, the frequency values tend to decrease for all three modes across all methods, and both L h 0 = 20 and 5.

Effect of rotating speed and hub radius
Table 3 presents the first three dimensionless natural frequencies of a rotating double-tapered homogeneous (g x = g z = 0) beam analyzed using three different theories: TSDT, FSDT, and CBT.The fre- quencies obtained from these theories are compared with the reference values from 23 using FSDT and 3D-FE results.The analysis is conducted with two different hub radii (R = 0 & L) , and three rotating speeds (� = 0, 50 and 100 rad/s) .The beam dimensions and material properties used in the analysis are 2 GPa, and ν = 0.3) .The ( 39) Table 1.Convergence of the first five dimensionless natural frequencies of a rotating BFGM double tapered beam.

Modes
Beam N e  3 reveal a notable level of agreement among the TSDT, FSDT, CBT, and 3D-FE methods with the reference results.Particularly, the CBT method exhibits excellent agreement in this specific example, where the beam under investigation is characterized by being extremely thin.

3D-FE
Regarding the effect of rotating speed, it is evident that an increase in rotational velocity leads to an increase in the natural frequencies.Furthermore, it can be inferred that the eigenfrequencies exhibit an increasing trend with the rise in angular velocity and hub radius.This behavior can be attributed to the centrifugal stiffness effect.As the angular velocity and hub radius increase, the centrifugal force acting on the beam also increases, resulting in a higher stiffness.Consequently, the eigenfrequencies of the beam are influenced and exhibit an upward trend in response.

New results
This subsection presents the vibrational behavior of a BFGM double tapered rotating beam.The beam material properties are E 0 = 210 GPa, ρ 0 = 7850kg/m 3 , and ν = 0.3 , and the beam dimensions are b 0 = 0.04 m and L = 0.4 m.Several parameters such as the index parameters ( g x and g z ) for the material distribution, rotating speed ( ), hub radius (R), and aspect ratio ( L h 0 ) are considered.

Effect of aspect ratio on natural frequencies: CBT, FSDT, and TSDT comparison
Figure 5 presents the variation of the first four dimensionless natural frequencies of a non-rotating (� = 0 rad/s) homogeneous g x = g z = 0 uniform (C b = C h = 0) beam versus the aspect ratio.The results demonstrate a correlation between the aspect ratio and the natural frequency, indicating that an increase in the aspect ratio leads to a higher natural frequency.However, this increase becomes less significant for the flapping modes as the aspect ratio reaches higher values.Furthermore, there is a deviation between the CBT theory and the other two theories, FSDT and TSDT, which is more pronounced for higher modes of flapping vibration.However, in the case of axial vibration, there is a good agreement between the three theories.Additionally, there is a good agreement between the FSDT and TSDT theories, particularly for low aspect ratios or thick beams.Moreover, with changing aspect ratio, the axial and flapping modes may interchange their order.This means that as the aspect ratio varies, there is a possibility for the axial mode and the flapping mode to switch positions in the mode sequence, resulting in a change in their relative order.

Effect of taper ratio on natural frequencies: CBT, FSDT, and TSDT comparison
Figure 6 presents the variation of the first four dimensionless natural frequencies of a non-rotating (� = 0 rad/s) homogeneous g x = g z = 0 beam versus taper ratios for aspect ratio L h 0 = 10 .Subfigures (a) and (d) depict the first flapping and axial modes, respectively.These subfigures provide clear evidence of a direct relationship between the taper ratio and the natural frequency, indicating that an increase in the taper ratio leads to an increase in the natural frequency.Furthermore, the highest natural frequencies are achieved when the taper ratio is increased in both the height and width directions, followed by tapering in the width direction only and then tapering in the height direction only.Notably, when the taper ratio is changed in either the height or width direction alone, the first axial mode exhibits identical natural frequencies due to the beam's consistent width and height.Moreover, a slight deviation can be observed in the results obtained from the CBT method compared to both the FSDT and TSDT methods in the first flapping mode.
In Subfigures (b) and (c), the second and third flapping modes are presented, respectively.These subfigures clearly demonstrate that an increase in the taper ratio in the width direction corresponds to an increase in the natural frequency.Conversely, increasing the taper ratio in both the height and width directions or in the height direction alone leads to a decrease in the natural frequency.Furthermore, the highest natural frequencies are obtained when the taper ratio is increased in the width direction only, followed by tapering in both the height and width directions and then tapering in the height direction only.Additionally, the observed deviations in the results obtained from the CBT method, compared to both the FSDT and TSDT methods, are more pronounced in the third mode than in the second mode.Table 3.The first three natural frequencies of a homogenous rotating double tapered for various rotating speed, and hub radii.7 presents the variation of the first four dimensionless natural frequencies of a homogeneous (g x = g z = 0) uniform (C b = C h = 0) beam versus rotating speed and hub radius for aspect ratio L h 0 = 10 .Figure 7a-c show the first three flapping modes, respectively.Increasing the rotating speed is associated with higher natural frequencies in all three flapping modes, attributed to the stiffening effect of centrifugal force.The hub radius also influences the natural frequency, particularly at higher rotating speeds, where the larger hub radius induces a stronger centrifugal force.Both increasing rotating speed and larger hub radius contribute to higher natural frequencies due to the enhanced effect of centrifugal force as shown in Eq. (24).Furthermore, the deviation between the CBT method and the FSDT and TSDT methods increases with higher modes, indicating that differences in predictions become more pronounced as the mode number increases.
Figure 7d indicates that the rotating speed and hub radius do not have a significant effect on the natural frequency in the first axial mode.

Effect of gradient index on natural frequencies: TSDT
The findings presented in Fig. 8 demonstrate the relationship between the gradient indexes g x and g z and the first four dimensionless natural frequencies of a uniform beam for two dimensionless speeds ¯ = 0 and 5 according to TSDT. Figure 8a-c results indicate a clear inverse correlation, where the dimensionless natural frequencies decrease as the gradient index increases.Furthermore, it is observed that the gradient index in the x-direction ( g x ) has a more significant impact on the natural frequencies compared to the gradient index in the z-direction ( g z ).The natural frequencies of all three flapping modes exhibit an increase as the rotating speed is increased.Figure 8d illustrates that the natural frequency in the first axial mode remains largely unaffected by changes in the rotating speed.

Uncertainties and stochastic results
Figures 9, 10, 11 and 12 show the Coefficient of Variation (C.O.V.) for the 1st and 2nd natural frequencies, denoted as V ω i (i = 1, 2) , of a uniform beam with R = 0.1 and L h 0 = 10 .These plots are generated as a function of rotational velocity.Figure 9 illustrates the Coefficient of Variation (C.O.V.) for the case where only the rotational velocity is treated as a random variable (Case 1), with a standard deviation of 10% .The analysis is conducted on a homogeneous beam ( g x = g z = 0 ) using the three different theories.Strong agreement is observed between the second-order perturbation and Monte Carlo methods, confirming the accuracy of the results across all rotational speeds.It is evident that V ω i increases with an increase in rotational velocity.Consequently, the C.O.V. of the first mode natural frequency remains nearly identical across all three theories, with only a slight deviation observed in the second mode, particularly between CBT and the other two theories, FSDT and TSDT.
Figure 10 represents the Coefficient of Variation (C.O.V.) for the case where only the material properties are considered as random variables (Case 2), with a standard deviation of 10% for E, G, and ρ .The analysis is per- formed on a homogeneous beam ( g x = g z = 0 ) using the three theories.It is observed that V ω i decreases with an increase in rotational velocity.Consequently, the C.O.V. of the natural frequency in the case of CBT deviates from that of FSDT and TSDT, with this deviation increasing from the first to the second mode.A comparison between Figures 9 and 10 reveals that the randomness of material properties has a more significant influence on V ω i than the randomness of rotational velocity within the speed range of 0 to 600 rad/s.
Figures 11 and 12 represent the Coefficient of Variation (C.O.V.) for the case where only the material distribution is subject to randomness (Case 3) according to TSDT. Figure 11a and b is applied to a BFGM beam ( g x = g z = [0.5 1 2] ) with a standard deviation of 10% for g x and g z .It is observed that V ω i decreases as the rotational velocity increases.Additionally, V ω i exhibits higher values for larger gradient indices.Figure 12a and b is applied to a BFGM beam ( g x = g z = 0.5 ) with standard deviations of [10 20 30 40 ]% for g x and g z .It is observed that V ω i decreases as the rotational velocity increases.Additionally, V ω i exhibits higher values for larger standard deviations.

Conclusions
In this study, a comparative analysis of free vibration behavior and stochastic analysis was conducted for rotating double-tapered beams composed of BFGM using three different beam theories: CBT, FSDT, TSDT.The material properties of the beams were characterized by an exponential distribution model.The investigation aimed to assess the impact of various parameters such as material distribution, taper ratios, aspect ratio, hub radius, and rotational speed on the natural frequencies.It also focused on studying uncertainties in material properties, material distribution, and rotational velocity.The results obtained from the study demonstrated the convergence and validation of the proposed model by comparing it with 3D-FE simulations conducted using ANSYS and previously published research.The convergence analysis indicated that the accuracy of the results reached a satisfactory level for practical purposes.Additionally, the stochastic analysis outcomes, derived using   the mean-centered second-order perturbation approach, underwent validation through comparison with the Monte Carlo method.Furthermore, the comparison between the beam theories showed that FSDT and TSDT exhibited a greater agreement with the 3D-FE simulations, particularly for thick beams, compared to CBT.This finding highlights the importance of considering shear deformation effects, which become more significant in thick beam structures.The results of the study showed that the natural frequencies of the rotating beams increased with increasing taper ratio in the width direction, aspect ratio, rotational speed, and hub radius.The results also showed that the natural frequencies of the beams decreased with increasing gradient index.In the context of random parameters, distinct trends emerge in the coefficient of variation (C.O.V.) concerning rotational velocity and various sources of uncertainty.Specifically, an increase in rotational velocity results in an increasing C.O.V. for random velocity, while random material properties and material distribution exhibit a decreasing C.O.V. trend with increasing velocity.Notably, for the first mode, a consistent C.O.V. is observed across all three theories for random rotational velocity.However, in the case of random material properties, discernible deviations are noted for CBT compared to FSDT and TSDT. https://doi.org/10.1038/s41598-023-44411-0

Figure 5 .
Figure 5.The first four dimensionless natural frequencies versus aspect ratio.

Figure 6 .
Figure 6.The first four dimensionless natural frequencies of a non-rotating homogeneous beam versus taper ratios.

Figure 7 .
Figure 7.The first four dimensionless natural frequencies of a homogeneous uniform beam versus rotating speed and hub radius.

Figure 8 .
Figure 8.Comparison of the first four dimensionless natural frequencies between non-rotating and rotating uniform beams versus the gradient index.

Figure 9 .
Figure9.C.O.V. of 1st and 2nd natural frequencies for random rotating velocity (Case 1) versus velocity for three beam theories using the second-order perturbation method compared to Monte Carlo (MC).

Figure 10 .
Figure 10.C.O.V. of 1st and 2nd natural frequencies for random material properties (Case 2) versus rotating velocity for three beam theories using the second-order perturbation method.